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Q^i Particle shape and aggregation have a strong influence on the spectral profiles of 

^ I infrared phonon bands of solid dust grains. Calculating these effects is difficult due 

' to the often extreme refractive index values in these bands. In this paper, we use 

c/3 . the Discrete Dipole Approximation (DDA) and the T-matrix method to compute 

^ ' the absorption band profiles for simple clusters of touching spherical grains. We 

invest reasonable amounts of computation time in order to reach high dipole grid 
^ ■ resolutions and take high multi-polar orders into account, respectively. The infrared 

phonon bands of three different refractory materials of astrophysical relevance are 
considered - Silicon Carbide (SiC), Wustite (FeO) and Silicon Dioxide (Si02). We 
demonstrate that even though these materials display a range of material properties 
and therefore different strengths of the surface resonances, a complete convergence 
is obtained with none of the approaches. For the DDA, we find a strong depen- 
dence of the calculated band profiles on the exact dipole distribution within the 
aggregates, especially in the vicinity of the contact points between their spherical 
constituents. By applying a recently developed method to separate the material op- 
tical constants from the geometrical parameters in the DDA approach, we are able 
to demonstrate that the most critical material properties are those where the real 
part of the refractive index is much smaller than unity. 
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Fig. 1. The clusters we study; a linear cluster of five particles lin5, a seven particle 
cluster frac7 and a eight particle cluster sc8. 

Grain growth by aggregation is an important process in dense cosmic envi- 
ronments as well as in the Earth's atmosphere. Besides influencing dynamic 
properties it also changes the absorption and scattering properties of the solid 
dust particles for electromagnetic radiation (e.g. [1,2]). This is especially true 
in spectral regions where resonant absorption occurs, such as the phonon bands 
in the infrared (IR). It is quite well known that shape and aggregation effects 
actually determine the band profiles of such absorption and emission bands, 
which hinders e.g. the identification of particulate materials by their IR bands. 
Despite many detailed investigations the full understanding of the influence of 
grain aggregation is still lacking, due to the complicated nature of the problem 
(e.g. [3,4]). 

We have investigated different simplified cluster shapes, as a means to compare 
how the numerical codes which are often used in astrophysical calculations 
manage to converge. The clusters consist of 5, 7 and 8 particles and are of 
significantly different shapes. For the material properties we have considered 
three astrophysically relevant materials differing in their optical properties; 
Silicon Carbide (SiC), Wustite (FeO) and Sihcon Dioxide (Si02). 

The paper is structured such that Sect. 1 presents a description of the cluster 
shapes. Sect. 2 provides a description of the different materials considered and 
their optical properties. In Sect. 3 the numerical methods used for determining 
the extinction of the clusters are described. Sections 4 and 5 describe our re- 
sults for the cluster-of-spheres method and the DDA calculations, respectively, 
and describe the convergence problems encountered by the different methods. 
Sect. 6 offers our conclusions. 



1 Cluster Structures 

We consider clusters of identical touching spherical particles arranged in three 
different geometries: fractal, cubic, and linear. For a high precision of the 
calculations, we restrict the number of particles per cluster to less than 10. 



2 




Fig. 2. Transmission Electron Micrograph of an /3-SiC cluster produced in the lab- 
oratory by laser-induced pyrolysis [6]. 

We have selected three main geometries, namely a "snowflake I'st -order pre- 
fractal" cluster (fractal dimension D = In7/ln3 = 1.77, [5]), where one sphere 
is surrounded by six others along the positive and negative cartesian axes, a 
cluster of eight spheres arranged as a cube and a linear chain of five spheres. All 
of the clusters (Fig. 1) consist of homogeneous spheres with radii R = 10 nm 
and embedded in vacuum (or air). 

An SiC cluster consisting of nano-particles produced in the laboratory by 
laser-induced pyrolysis [6], which has served as inspiration for this paper, is 
shown in Fig. 1. 



2 Infrared properties of our materials 

In order to study the convergence of the different codes, we chose three dif- 
ferent materials of which our clusters are composed: crystalline /3-SiC, crys- 
talline FeO and amorphous Si02. These solids differ significantly from each 
other in their infrared properties. First, while /3-SiC and FeO have only one 
mid-infrared resonance, amorphous Si02 has three of them. Second - and more 
importantly - the minima and peak values of these materials' optical constants 
are very different from each other, see Fig. 2. Roughly speaking, the maximum 
n and k values of FeO, by reaching 5-6 close to 30 /im, amount to only 50% of 
those of /?-SiC in its resonance at 12.6 /xm, and the same holds true for Si02 
as compared to FeO. Due to its sharp resonance, /3-SiC not only reaches a 
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Fig. 3. The optical constants n and k for the three different materials considered for 
the clusters; SiC, FeO and Si02. It can be seen that SiC possesses a significantly 
stronger resonances than the two other materials. 

very high maximum, but also a very deep minimum of n, which makes the 
calculation of its absorption and scattering behaviour especially difficult (see 
also Sect. 4 and 5). More details on the maxima of n and k as well as the 
minima of n for SiC, FeO and Si02 can be found in Tab. 1. 

Table 1 

Maxima of the optical constants n and k of SiC, FeO and Si02. Amorphous Si02 
has three resonances in the mid-IR range which have been labeled rl . . . r3. 



Resonance 


^max 


^(j^max ) 




■^(llmin) 










[fim] 




[fim] 




[fim] 


/3-SiC 


13.0 


12.6 


0.12 


10.8 


12.6 


12.5 


FeO 


6.3 


33.0 


0.99 


20.3 


5.0 


29.2 


SiOa, rl 


2.8 


9.8 


0.36 


8.9 


2.6 


9.2 


SiOa, r2 


1.8 


12.8 


1.62 


11.9 


0.3 


12.4 


SiOa, r3 


2.9 


22.7 


0.58 


20.4 


2.4 


21.3 



The optical constants used for the amorphous Si02 we got from [7] and for FeO 
from [8]. The data for /9-SiC were calculated from a Lorentz oscillator model 



4 



A 
o 
V 



(U 

O 



5 



10 - 



jS-SiC: lin 5 cluster 
SCSMTM results 



doshed: L=7 




dotted: 
L=29 



0.5 



12.0 



12.5 



11.0 11.5 
X [/xm] 

Fig. 4. The calculated extinction of a SiC linear five particle cluster lin5 using the 
T-matrix SCSMTM code by [10]. The polar order, L, of each calculation is indicated 
by an annotation. 



according to [9], assuming a damping constant, 7, of 10cm~^. It should be 
noted that high-quality (terrestrial) SiC is characterized by damping constants 
of only 1-3 cm~^, implying an even sharper resonance than for 7 = 10cm~^. 
However, circumstellar SiC is almost certainly characterized by an imperfect 
lattice structure, corresponding to a damping constant significantly larger than 
3cm~^ (see [9] for more details). 



SiC, FeO and Si02 are not only among the potential condensates in circum- 
stellar shells, but due to the above indicated values of their infrared optical 
constants also well-suited for testing the convergence of the different codes 
delivering the cluster absorption and scattering effiencies. As we are going to 
show, covergence in the absorption efficiency calculations cannot be reached 
by means of the codes we used for materials with optical constants reaching as 
deep minima and as high maxima of n as /?-SiC does; but also for a material 
like amorphous Si02 the convergence of our calculations is very slow. 
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Fig. 5. The calculated extinction of a FeO linear five particles cluster lin5 using the 
SCSMTM code [10]. The polar orders shown span from L is 7 to 31. 



3 Computational approaches 



The clusters-of-spheres calculations have been performed using (1) the T- 
matrix code by D.W. Mackowski (SCSMTM) calculating the random-orienta- 
tion scattering matrix for an ensemble of spheres [10] and (2) the program 
developed by M. Quinten [11] (MQAGGR, commercially available) based on 
the theoretical approach by [12]. Both programs aim at solving the scattering 
problem in an exact way by treating the superposition of incident and all scat- 
tered fields, developed into a series of vector spherical harmonics. Available 
computer power, however, forces to truncate the series at a certain maximum 
multipolar order L = npolmax, which in both programs can be specified ex- 
plicitly. Furthermore, both programs perform an orientation average of the 
cluster, the resolution of which was set to 15 degrees in MQAGGR and 10 
degrees in SCSMTM for theta (the scattering angle). The variation in the azi- 
muthal angle is not specified in SCSMTM. In MQAGGR it is varied between 
and 360 degrees, again with a resolution of 15 degrees. 

The discrete dipole approximation (DDA) method is one of several discretiza- 
tion methods (e.g. [13,14]) for solving scattering problems in the presence of a 
target with arbitrary geometry. In this work we use the DDSCAT code version 
6.1 by B.T.Draine and P.J. Flat au [15], which is very popular among astro- 
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Fig. 6. The calculated extinction of a Si02 linear five particles cluster lin5 using the 
SCSMTM code [10]. The polar order L shown spans from 7 to 31. 



physicists. In DDSCAT the considered grain/cluster is replaced by a cubic 
array of point dipoles of certain polarizabilities [16]. The cubic array has nu- 
merical advantages because it allows a significant speedup of the conjugate 
gradient method, applied to solve the matrix equations describing the dipole 
interactions, by using Fast Fourier Transforms. By specifying an appropriate 
grid resolution, calculations of the scattering and absorption of light by inho- 
mogeneous media such as particle aggregates can be carried out to in principle 
whatever accuracy is required. The size of the dipole grid is in reality limited 
by the available computing power. 

M. Min et al. [17] have presented an easy to use method to compute the absorp- 
tion and scattering properties of small particles with arbitrary shape, struc- 
ture, orientation and composition based on a solution of the DDA equations 
in the Rayleigh domain. For a given geometrical shape of the particles, the 
solution has to be computed only once to obtain the absorption and scatter- 
ing properties for arbitrary values of the refractive index. This method thus 
provides a nice tool to study the dependence of the convergence behavior of 
the DDA on the value of the refractive index of the particle. 
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4 Convergence of the clusters-of-spheres method 

For a pair of gold spheres, it has been shown already by [4] that the clusters-of- 
spheres method, when applied to touching particles, hardly yields convergent 
results when including multipolar terms up to an order of eight. By comparing 
their results - obtained according to the formalism by [12] - with the exact 
quasistatic solution for a pair of touching particles. Smith et al. found that 
especially for the longitudinal optical mode, the multipolar expansion up to L 
= 8 does not lead to satisfactory results. (Note that the values of n and k for 
gold near the 484 nm plasmon resonance peak studied by [4] amount to n^es 
= 1.1 and kres — 1-8. Beyond Xres, k(A) increases with A.) 

In our application of the clusters-of-spheres method, we came to similar con- 
clusions, even though we included multipolar terms up to an order of L = 31. 
This is illustrated in Figs. 4-6, which show the extinction band profiles for 
the linS cluster of SiC, FeO and Si02 spheres in the 10-12.5, 15-32 and 8.3- 
9.5 /xm wavelength ranges, respectively. The spectra are presented in terms of 
the extinction efficiency Q,ext=Cext/i'^<^>'^), with C^xt being the extinction 
cross section and <a> the radius of a sphere with the same volume as the 
cluster. Qext is normalized by <a> to become independent of <a> within 
the Rayleigh limit. In the figures we plot Q,ext/<^> because this quantity is 
independent of the particle size when the particles are in the Rayleigh limit. 

For all three materials, a phonon band - split up into at least two components 
- is seen in the extinction spectrum. Also in all three cases, for the 'blue' 
component of the phonon band, convergence is achieved for L > 20, while for 
the 'red' component(s), this is not the case, but instead, even for L = 31, there 
is a continuous peak shift if L is increased. Recall that our extinction spectra 
do not refer to specific orientations of the particle chain with respect to the 
external electromagnetic field, but are orientationally averaged. By studying 
specific orientations of the chains, it becomes evident that the 'blue' phonon 
band component corresponds to an orientation of the electric field perpendicu- 
lar to the long axis of the chain, while the 'red' component corresponds to the 
complementary case. It is this latter case which causes the slow convergence. 
The convergence, however, also depends on the values of optical constants. 
For amorphous Si02, with its smaller maximum values of n and k (cf. Tab. 1), 
the differences between the results for the individual multipolar orders are 
significantly smaller than for FeO with its larger peak values of n and k. 

The 20 /im phonon band of amorphous Si02 (not shown in our figures) actually 

converges even faster than the the 9 fim band. Here, convergence seems to be 
achieved for L<31. This is probably due to the fact that Re(m) does not reach 
as low values as for the 9 /^m band. This dependence of the convergence on 
Re(m) is also observed with the DDA approach and will further be discussed 
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Table 2 



The number of dipoles used with the DDSCAT [15] calculation as well as the number 
of dipoles actually within the cluster. 



Cluster 


# dipoles 


# dipoles 
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# dipoles 


# dipoles 




in the grid 


in the cluster 




in the grid 


in the cluster 


ill iU 
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94y94y94 — 1 "^894 


71 60 


lin5 


60x60x60 = 216000 


4475 


sc8 


25x25x25 = 15625 


8217 


lin5 


72x48x32 = 110592 


7792 


sc8 


60x60x60 = 216000 


111976 


lin5 


80x80x80 = 512000 


10515 


frac7 


24x24x24 = 13824 


1757 


lin5 


90x60x40 = 216000 


14845 


fracl 


25x25x25 = 15625 


2105 


lin5 
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18620 


fracl 


60x60x60 = 216000 


28973 




Fig. 7. The calculated extinction of a SiC linear five particle cluster lin5 using the 
DDSCAT code by [15] . The number of the dipoles within the cluster is indicated by 
the annotations, see also Tab. 2. For comparison the result as calculated with the 
SCSMTM code [10] for L = 31 is shown. 



in the next section. 
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Fig. 8. The calculated extinction of a FeO linear five particle cluster lin5 using the 
DDSCAT code [15]. The number of the dipoles within the cluster is indicated by 
the annotations, see also Tab. 2. For comparison the result as calculated with the 
SCSMTM code [10] for L = 31 is shown. 

5 Convergence of the DDA methods 

As well as the cluster-of-spheres methods, the DDA method also encounters 
difficulties with reaching convergence for the clusters we have chosen to study, 
see Figs. 7 - 10. For the /m5 cluster composed of SiC it is clear from Fig. 7 that 
we are most likely not even close to obtaining a converged result, although 
we have used rather large numbers of dipoles (see Tab. 2) for the calculations. 
This is especially true at wavelengths longward of the primary resonance at 
10.6 fim. Some of the spectra obtained show secondary resonances similar to 
the SCSMTM results whereas others are relatively flat. Similarly, for the FeO 
lin5 cluster the results seem to be converged towards the short-wavelength 
part of the spectrum shown in Fig. 8, but not towards the longer wavelengths, 
especially for one particular spectrum (4475 dipoles) which strongly differs 
from all others including the SCSMTM resuh. For SiOa (Fig. 9) the DDSCAT 
calculation is closer to converging due to the lower values of the optical con- 
stants for this material as compared to FeO and SiC, see Tab. 1. However, 
again the 4475 dipole spectrum shows a secondary resonance which is neither 
present in the other three nor in the SCSMTM result. Similar results concern- 
ing the convergence were found for the two other cluster shapes studied frac7 
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Fig. 9. The calculated extinction of a Si02 linear five particle cluster lin5 using the 
DDSCAT code [15]. The number of the dipoles within the cluster is indicated by 
the annotations, see also Tab. 2. For comparison the result as calculated with the 
SCSMTM code [10] for L = 31 is shown. 

and sc8, Fig. 10. 

An influence of the refractive index values on the precision of the DDA is 
expected. For materials with large refractive index (|m| > 2) [18] show that 
especially the absorption is overestimated by DDA. According to [18], the 
lattice dispersion relation prescription for the polarisability, ««, gives fair ac- 
curacy for scattering but poorer results for absorption. When |m — 1| is large, 
the continuum material is effective at screening the external field: in the limit 
|m — 1| — > oo the internal field generated by the polarization would exactly 
cancel the incident field, so that the continuum material in the interior of the 
target would be subjected to zero field. In the case of a discrete dipole array, 
the dipoles in the interior will also be effectively shielded, while the dipoles 
located on the target surface are not fully shielded and, as a result, absorb 
energy from the external field at an excessive rate. This error can be reduced 
to any desired level by increasing the number of dipoles, thereby minimizing 
the fraction N~^^^ of the dipoles which are at surface sites, but very large 
values of N are required when |m — 1| is large [18]. 

At a closer view to the results of the DDA calculations, we found that the 
choice of the dipole grid is significant for the obtained result. For the linear 
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Fig. 10. The calculated extinction of the SiC clusters frac7 and sc8 using the 
DDSCAT code [15]. The number of dipoles within the cluster is indicated by the 
annotation, see Tab. 2 for details. 

five particle cluster lin5, a grid with a dipole number which is a multiple of 
five provides a significantly different result than dipole grids which are not 
a multiple of five. The former are the ones which generally produce stronger 
secondary resonances as can be seen from a comparison of the resulting spectra 
(shown in Fig. 7) with Tab. 2 (see also Fig. 11, which will be discussed later). 
For the semi-fractal seven particle cluster frac7 the difference was between 
dipole grids which were a multiple of three and for the cubic eight cluster sc8 
it was for dipole grids which were a multiple of two, equivalent to the length 
scale of the longest chain within the cluster. 

The results shown in Figs. 7-10 were calculated using DDSCAT [15] but 
identical results were found using the code by [17]. According to this solution 
the absorption cross section averaged over all orientations of an arbitrarily 
shaped particle obtained using DBA can be written as 

3^ Wj 

C'abs — 



kV Im ( ) 



where k = 271 / X, V is the material volume of the particle, is the number of 
dipoles. The > > 1, the so-called form-factors (not to be confused with 
the symbol for the multipolar order), and Wj, the weights, are obtained from 
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Fig. 11. A demonstration of how the results of the DDA calculations for the linS 
cluster show a dependency on the choice of dipole grid used. The top panel shows 
the SiC cluster, the middle panel the FeO cluster and the bottom panel the Si02 
cluster. The results are obtained with the code by [17], an identical dependence is 
found by the DDSCAT code of [15]. The dark gray area show the range of seven 
spectra obtained when using a dipole grid which is a multiple of five. The light gray 
shaded area is the range of 24 spectra obtained when using dipole grid which isn't 
a multiple of five. A similar result is found for the frac7 and sc8 clusters but here 
for grids which are multiple of three and two, respectively. 

the DDA equations. Once the Lj and the Wj are found, it is trivial to compute 
the absorption cross section for any wavelength or complex refractive index, 
m. 

Fig. 11 shows a selection of results for the linS cluster, divided in into grids 
dividable and non-dividable by 5. It is obvious, especially for SiC and FeO, 

that the dividable grids produce secondary resonances, which lead to long- 
wavelength extinction enhancement but with a large spread. Consequently, 
they are the ones that contribute stronger to the convergence problems than 
the non-dividable grids. 
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Fig. 12. Contour plots of Cabs/(^^) as a function of the real and imaginary part 
of the refractive index for the linS cluster, averaged over all particle orientations. 
The contours are obtained using DDA with a dipole grid which has 50 (left panels) , 
52 and 54 (middle panels) and 55 (right panels) dipoles along the long axis of the 
system. The lower panels arc simply blowups of the lower left corners of the upper 
panels. The contours are plotted with increments of 0.5. 

The method by [17] allows to examine these effects in even more detail by 
investigating the origin of the resonances. According to Eq. 1, each form- factor 
Lj causes a resonance in the complex m plane at 

™=^ri. (2) 



In the limiting case ^ oo the Lj and Wj form a continuous distribution, 
smoothing out all these resonances to a continuum function of m. However, 
using a limited number of dipoles will result in a limited number of resonances. 
Especially for small values of Re(m) , this resonance behavior might show up. 
The positions of these resonances are sensitive to the exact configuration of 
the dipoles in the particle. Thus it is very hard to obtain convergence using 
DDA for low values of Re(m). 

The previous is illustrated in Fig. 12 where we plot in the contours the dimen- 
sionless quantity Cabs/ {kV) for the lin5 aggregate as a function of the real and 
imaginary part of the refractive index. We do not give any numerical values 
for the contour lines because we would like to focus on the general shape of 
the contours. In the upper panels it can be seen that there is a large differ- 
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ence between the contours obtained using different numbers of dipoles. It is 
clear in all plots that a number of resonances are present along the axis with 
Re(m) = 0. More generally speaking, the contours look very different over the 
whole complex m space. If we zoom in on the region with moderate values of m 
(the lower panels), we see that the differences between the 52 and the 54 wide 
dipole grids are almost gone in this region of the complex m plane. All the 
differences between the contours shown in the middle two panels are located 
at very small values of Re(m). However, the grid of 50 and 55 dipoles width 
still displays a different shape of the contours. More generally, for this cluster 
shape there is a different behaviour when the number of dipoles along the long 
axis of the system is a multiple of five. This illustrates the bad convergence 
behavior of this cluster shape when using DDA. 



6 Conclusions 

We have studied the performance of two cluster-of-spheres ([10,11]) and two 
discrete dipole approximation (DDA; [15,17]) methods for calculating the ex- 
tinction of aggregates. The methods were chosen because they are popular and 
often used in studies of astrophysical dust extinction and scattering problems. 

We present results of the calculated extinction within infrared absorption 
bands of SiC, FeO and Si02 clusters composed of 5 to 8 spherical particles 
of radii 10 nm. The clusters had three different shapes: hnear chain, semi- 
fractal and simple cubic. The materials display a range of material properties 
and therefore have different strengths of the surface resonances. When the 
real part of the refractive index is much smaller than unity, none of the four 
methods are able to converge. For the two DDA methods there is a strong 
dependence of the calculated band profiles on the exact dipole distribution 
within the particles. For the linear five particle cluster lin5 the result depend 
on whether the grid is a multiple of five or not, for the semi-fractal seven parti- 
cle cluster frac7 and the simple cubic eight particle cluster .sc5thc dependence 
is on whether the grid is a multiple of three and two, respectively. 

Currently it does not seem possible to calculate the absorption efficiencies of 
clusters of spheres for materials with optical constants n << 1 in a reliable 
way. We assume that the critical point is the contact geometry of the par- 
ticles, which explains the importance of the resolution of the methods used. 
This implies that for corresponding experiments it is of great importance to 
study the contact points between the particles. Without a good investigation 
of the contact points a comparison between experimental data and model cal- 
culations are very likely to be fruitless. One could hope that things might be 
easier when particles are less perfect than the spheres we have considered, 
since that will give rise to fewer resonance effects. 
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